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Abstract 

A finite sum of exponential functions may be expressed by a linear 
combination of powers of the independent variable and by successive in- 
tegrals of the sum. This is proved for the general case and the connection 
between the parameters in the sum and the coefficients in the linear com- 
bination is highlighted. The fitting of exponential functions to a given 
data- set is therefore reduced to a multilinear approximation procedure. 
The results of this approximation do not only provide the necessary in- 
formation to compute the factors in the exponents and the weights of the 
exponential terms but also they are used to estimate the errors in the 
factors. 



1 Introduction 

From time to time the need arises to fit a sum of exponentials to numerical 
data. That means to approximate a given data- set consisting of pairs of real 
numbers (Xj,yj) by the following expression: 

N 

yixj^ao + ^a^* (1) 

i=l 

where x, Xj G M + , yj G M, N G N + and a;, bi are unknown real numbers which 
have to be chosen so that the fit becomes optimal. 

If the bi were known, the task usually would be a well posed linear problem, 
but if the bi are unknown too, it turns out to be ill conditioned. The hopelessness 
of efforts dealing with this kind of problem has been described drastically by 
F.S. Acton ^] in a chapter entitled "What not to compute". 

At first sight, fitting equation ^ to a given data- set inevitably seems to be 
a nonlinear problem. However it has been noted [3] that equation ^ may 
be expressed as a linear combination of powers of x and successive integrals of 
y(x), reducing the problem to a multilinear fitting procedure. This method is 
based on the fact that y(x) can be shown to satisfy an ordinary linear differential 
equation of N-th order with constant coefficients. The roots of the characteristic 

'Metallurgical research department, voestalpine, 4031 Linz, Austria, 
e-mail: bernhard.kaufmann@voestalpine.com 



1 



polynomial of this equation give the bi and the at are identified as solutions to 
linear equations involving the bi and the derivatives of y(x) at x = 0. However, 
derivatives of experimental data- sets enhance the errors in the data, therefore 
it is desirable to eliminate them. Actually, [3] shows already the way to do this, 
but only for the case N = 2 and ao = 0, the general case not really being obvious. 
In the present paper, this method of linearizing the fitting procedure is revived 
and without referring to differential equations and without using derivatives, 
the general case is derived. Additionally, a method for estimating the errors in 
the exponential factors is presented. Of course, the problem remains ill posed, 
but linear fitting offers computational advantages over nonlinear approximation 
and also supplies estimates of the errors in the computed coefficients, which may 
be used to predict the errors in the exponential coefficients. 



2 Results 

The functions used to construct a linear approximation problem are powers of x 
and successive integrals of y(x). Before the announced relation can be asserted, 
some definitions are required. 

Definition: The k-th integral of y(x) is defined recursively: 

= y(x) 

h(x) = [ h-i(t)dt k>0 (2) 



Definition: ftNij, otNij- Given the set B = {bi, • ■ • , 6jv} of N exponential 
factors bk in equation consider the products of i different elements of B. 
Each of these products corresponds to a combination of i elements out of B, the 
number of these products therefore is 

AT! / N\ 

as is proved in combinatorics. We assume that the products are ordered in some 
way. Pwij then is the j-th of these products. Additionally, we define /3nqi = 1. 

cxNij is the sum of all ak in equation ^ excluding those whose index is equal 
to that of one of the b's in Pmj- By definition, 01^01 = X)fclo afe - Obviously, 
each ctNij contains at least ao. 
Example: N=3 



2 



/?301 


= 1 










/Sail 


= 6l 


0312 


= b 2 


0313 


= b 3 


/?321 


= &1&2 


/?322 


= bib 3 


@323 


= &2&3 


/?331 


= bib 2 b 3 











"301 = a + ai + a 2 + a 3 

"311 = a + a 2 + a 3 a 3i2 = a + ai + a 3 a 313 = a + ai + a 2 

"321 =0-0 + 0,3 a 322 = a + a 2 a 323 = a + 01 

"331 = a o 

With these definitions, the central statement of this article now may be 
asserted: 



jV N C(N.i) N l C(N,i) 

a o + ^a t e- b > x = -^2,Ii(x) Pnij+Y— Y I^NijOiNij (4) 

i— 1 i— 1 j— 1 i—0 j — 1 

Assuming the validity of equation 0] the task now consists in approximating 
the data- set {(xj,yj)} by a linear combination of the 2N functions (Ii, . . . ,I N , 
x, . . . ,x N ) plus a constant. By standard linear approximation techniques the 
coefficients (ci, . . . , cat, d\, . . . , d/v) and the intercept do may be determined 
together with their errors (Aci, . . . , Acjv, Adi, . . . , AcIn) and Ado- It follows 
that 



C(N,i) 

d = - Pmj i = l,...,N (5) 

C(N,i) 

di = Y V ' y ,' V " i = l,...,JV (6) 

iV 

do = "Jvoi = 2J ak ( 7 ) 

fc=0 

Given the Cj inElVieta's root theorem asserts that the 6$ are the N roots of 
the polynomial 

N 

P(x) = x N + J2(- 1 Y +1 ^ N ~ t (8) 

i=l 

As soon as the bi are known, the expressions El and [7| represent a system of N+l 
linear equations for the N+l coefficients tij. 

If the Aci are small, the relation between the errors may be approximated 
by the linear terms of the Taylor- series for P(x). 

ATM + 

2—1 
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As the bk are roots of P, AP should be zero and therefore, inserting bk for 
x, we get: 



/V 



dP(b k ) E 
dx i=l 



dP(b k ) 



Ac- 



tio) 



Treating the Cj and bk as probability variables with standard deviations s C i and 
Sbk, the standard deviation and therefore the estimated error of bk is given by 



S&fc 



1 



dx 



N 



N N 



(ii) 



As usual, Cov(a, Cj) means the covariance between a and Cj. 

It remains to show that equation 01 is valid. For this purpose it is useful to 
state some properties of the coefficients (3. 

For any I with 1 < I < N and any i < N the sum of all Pniiu may be divided 
into the sum of all /3/v; TO containing bi and those not containing bf. 



C(N,l) 

E 

r. 



C{N-l,l-l) 



C(JV-1,;) 



E P(N-l)(l-l)j + E Pff-1,1,, 



(12) 



m — 1 



3=1 



With Pny— 1)0 we denote the products not containing bi that is, which are 

chosen from the set {&1, . . . , £>i_i, &i+i, ■ ■ • , 6 at} containing N-l elements and 

f— i) 

not containing 6^. For I < we define Pny-iMm — 1- 

An important special case ofEl resu lts if i = N. Then Wjy-iMm = fi(N-i)im 
and the following expression results: 



C(AT,0 C(N-l,l) C(N-l.l-l) 

E 0Nljm= E 0{N-l)lj + bN E P{N-l){l-l)j (13) 



m— 1 



For a proof of equation 0] consider the following system of equations: 



4(z) 



AT 

a + E a * e 

«o ¥ + E 

i— 1 



-biX 



(14) 



The validity of 1141 is easily seen by performing the integrals in equation [21 
analytically. 

Now consider the following linear transformations defined recursively on the 
set of equations 1141 
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4 1} = h + hh +1 

t ] = it^ + bulfr^ h>l (15) 
For this kind of transformation a rather general relationship holds: 

h C(h,l) 

(=1 m=l 

Proof: Induction for h. For h = 1, propositionllfiliust repeats the definition 
of4 X) . Now assume that US] holds for . Then 

T (h+i) _ T (h) , Ah) _ 

h C(h,l) h C(h,l) 

h + ^2 ^ Phlmlk+l + bh+llk+l + ^ ^ bh+lf3hlmlk+l + l = 
l — l m— 1 l — l m— 1 

h C(h.l) h + lC(h,l~l) 

Ik+^^2 Phlmh+l + bh+llk+l + ^ bh+l^h(l-l)mh+l = 

Z=l m=l /=2 m=l 

C(/»+l,l) 'i C(/»+l,0 

4+ ^ /3(fe+i)im-^fe+i + ^ ^ (3(h+i)imh+i + bh+iPhhih+h+i 

m— 1 /— 2 m— 1 

where IT31 has been used in order to obtain the last line. Obviously, this result 
may be converted into 

h+lC(h+l,l) 

= J k + ^2 ^2 P{h+\)lmh+l 
1=1 m=l 

whereby the proof of El is completed. 
Consider now Jq \ 

JV C(N,l) 

f =y + EE feJi (17) 

;=i m=l 

Inserting 1141 this expands into 

iv c(jv,o 



y^ y^ Pmmii — 



N C(N,l) 

E E /w 

Z=l m=l 



iV i-l 



'"7T E"-E' 6 *> 



(18) 

The motive for applying transform ^| to I was to get rid of the exponential 
terms. The following proposition asserts that equation ED is actually free of 
exponential terms: 



h C(h,l) , y 

E E ^™( — ) = _1 l <i<h (19) 

1 = 1 m—l \ 1 / 

Proof: 

For h = 1 the assertion is trivial. Now assume that El is valid for h. Then 
the following calculations prove the truth for h+1 and therefore for all h: 



h+l C(h+l,l) 

E E P{h+l)lm 
1=1 m=l 



1 



E^ 



C(h,i) C(h,J-l) 

E Phlm+bh+l E Ph(l-l)m 
m—l m—l 



-1 



h+l 



h+l)l 



h h-lC(hA) t+1 

-i-^+^EE^fv) + 



&ft+i ^h+i 



-1 m—l 

h C{h,l) 



h+l 



P(h+l)(h+l)l 



EE a 



hint 



1=1 m=l 



-iV , b h+1 ( v 



ft.+i 

/3(/i+i)(/i+i)i 

Using equation El for and collecting all terms the last expression evaluates to 
-1. 

To complete the proof of equation 0] some more transformations on formula 
El are required: 



N C(N,l) N C(N,l) 

V + E E PNlmh = «0 + E E @ Nlr 



1 = 1 m=l 



1=1 m=l 



I N i-l j 

a °fr-E ai E(-^) J ^7r 



N N 

E a *+E 

i=0 1=1 



C(N,l) N C(N,j) N 

E 0Nlm,aQ - &Njm E ^(^i) 1 ' 3 

m—l m—l i— 1 



G 



N N 

E a *+E 

i=0 1=1 



C(N,l) N I C(N,j) N 

E PNimao + ^2ai(~bi) 1 E PNjm^CLii-bi) 1 ^ 

i=l j=l 7n— 1 2= 1 



m— 1 



w Z 
X 



i=0 Z=l 



E f^Nhn E fl « + E ai ^ iZ 
m=l i=0 i=l 



(20) 



where 



i-i C(N,l-p) 

5 "=E E ^(l-p)m(-6i) p + for Z > 1 

p=l m=l 



and 



5a = —6,: for Z = 1 
Using El for Z > 1 transforms into 

Su = 



l-l (C{N-l,l-p) 

E E a 

p— 1 \ 771=1 



(-;) 

(JV-l)(/-p)m 



C(iV-l,i-p-l) 



Substituting in the second part of this sum q for p+1 this expression transforms 
into 



C{N~l,l-l) 



l-l C(N-1A- 



=2 m=l 



S*l — E ^( AT-lia-1 wA + E E ^(JV-l)(i-p)m( 

m=l 

j_l C(N-l,l- 

-E E 

9 =2 m=l 



E E /V-i)(<- 9)m (-k) 9 - (-*)' + (^y 



C(N-X,l-1) 



E a 



C-i) 



(N — l)(l — l)m 1 



m=l 



Therefore Si; is the negative sum of all /Jrvz-m which contain hi. Consequently, in 
expression 1201 only those a% are not cancelled for which /3jv; m does not contain 
For that, [2U] may be written as 

N C{N,l) N j C(N,l) 

V + E E P Nlm ^ 1 = E 7T E PNlmOtNlm 



1 = 1 m=l 



Z! 

1=0 m=l 



which proves equation 0] 

Example: Consider this sum of two exponentials and a constant: 



y{x) = 0.3 + exp(-0.7x) + 0Aexp(-0.3x) 
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The function is evaluated in in the interval < x < 6 at Np equally spaced 
points. The discrete function values are multiplied by one plus a gaussian dis- 
tributed random variable so that the relative error has the standard deviation 
a. II and 12 are calculated using the trapezoidal method. For different settings 
of Np and a the coefficients cl and c2, dO (the intercept), dl and d2 and the 
corresponding errors of cl and c2 as well as the covariance between these two 
factors are determined by the commercial statistics program STATISTICA® 
and are listed in table 1. The parameters bl, b2, aO, al, a2 and the errors A&i 
and A6 2 are calculated from these coefficients as described above and are listed 
in table 2. 



Table 1: Statistically determined coefficients 



Id 


Np 


a 


Cl 


Aci 


C2 


Ac 2 


Cov 


do 


di 


d 2 


1 


601 


0.0000 


-1.0000 


io- b 


-0.2100 


10"° 





1.7000 


0.8800 


0.0315 


2 


601 


0.0001 


-1.0054 


0.0024 


-0.2130 


0.0014 


3.10~ 6 


1.7001 


0.8890 


0.0320 


3 


601 


0.001 


-1.0358 


0.0253 


-0.2300 


0.0146 


0.000343 


1.7008 


0.9396 


0.0349 


4 


601 


0.01 


-1.0842 


0.2954 


-0.2615 


0.1700 


0.0502 


1.7006 


1.0249 


0.0408 


5 


2001 


0.01 


-0.9155 


0.1190 


-0.1603 


0.0685 


0.00815 


1.6999 


0.7369 


0.0228 



Table 2: Parameter estimates based on the coefficients in table 1 



Id 


61 


A&i 


b 2 


Ab 2 


a 




0.2 


1 


0.7000 


3. 10-" 


0.3000 


3.10-° 


0.30 


1.00 


0.40 


2 


0.7020 


0.0019 


0.3034 


0.0021 


0.30 


0.99 


0.41 


3 


0.7134 


0.0180 


0.3224 


0.0196 


0.30 


0.95 


0.45 


4 


0.7220 


0.1211 


0.3622 


0.1754 


0.31 


0.88 


0.51 


5 


0.6796 


0.0281 


0.2359 


0.0911 


0.28 


1.09 


0.32 



The results show that for small errors in the coefficients the estimated vari- 
ance of bi and b 2 is also small and the estimate is realistic. The first case was 
computed without artificial noise, in this case the accuracy seems to be deter- 
mined mainly by the statistics program. Adding noise deteriorates the accuracy 
of the results rapidly. While a relative error of 0.0001 (case 2) still leads to 
a reasonable result, the tenfold relative error (case 3) already means that the 
calculated uncertainty of 6 2 is about 7%. A one- percent inaccuracy in the data 
(case 4) gives a result even with the first digit uncertain. As case 5 where the 
number of data- points is raised to 2001 shows, increasing the size of the data- 
set may at least partially compensate for noise. 
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